In [1]:
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline as IUS
from StellarSpectra.grid_tools import HDF5Interface
myInterface = HDF5Interface("../libraries/PHOENIX_submaster.hdf5")
In [2]:
def load_pixel(index, temp, logg, Z):
'''
Return the flux value at a specific pixel in a spectrum, defined by params
'''
params = {"temp":temp, "logg":logg, "Z":Z, "alpha":0.0}
flux = myInterface.load_flux(params)
return flux[index]
def get_wl(index):
return myInterface.wl[index]
In [4]:
#Determine truncation indices for order of interest
wl = myInterface.wl
inds = np.argwhere((wl > 5120) & (wl < 5220))[np.array([0, -1])]
inds[0], inds[1]
Out[4]:
In [37]:
flux = myInterface.load_flux({"temp":6000, "logg":4.0, "Z":0.0, "alpha":0.0})
print(len(flux))
#plt.plot(flux[inds[0]:inds[1]])
#plt.show()
In [6]:
temps = myInterface.points["temp"]
In [14]:
#Pixel indices for lines
cont = 5018
mg_bot = 5268
mg_side = 5242
fe_bot = 5053
fe_side = 5050
In [9]:
flux_cont = [load_pixel(cont, temp, 4.0, 0.0) for temp in temps]
flux_mg_bot = [load_pixel(mg_bot, temp, 4.0, 0.0) for temp in temps]
flux_mg_side = [load_pixel(mg_side, temp, 4.0, 0.0) for temp in temps]
flux_fe_bot = [load_pixel(fe_bot, temp, 4.0, 0.0) for temp in temps]
flux_fe_side = [load_pixel(fe_side, temp, 4.0, 0.0) for temp in temps]
In [11]:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.plot(temps, flux_cont, label="continuum")
ax.plot(temps, flux_cont, "ko")
ax.plot(temps, flux_mg_bot, label="Mg bottom")
ax.plot(temps, flux_mg_bot, "ko")
ax.plot(temps, flux_mg_side, label="Mg side")
ax.plot(temps, flux_mg_side, "ko")
ax.plot(temps, flux_fe_bot, label="Fe bottom")
ax.plot(temps, flux_fe_bot, "ko")
ax.plot(temps, flux_fe_side, label="Fe side")
ax.plot(temps, flux_fe_side, "ko")
#try fitting a spline to fe_side
myspline = IUS(temps, flux_fe_side)
fine_temps = np.linspace(5000, 7000, num=300)
fine_fe_side = myspline(fine_temps)
ax.plot(fine_temps, fine_fe_side, "k", lw=0.5)
ax.set_xlabel(r"Temperature ($K$)")
ax.set_ylabel(r"$\propto f_\lambda$")
ax.legend(loc="upper left")
fig.savefig("../plots/interpolation_temp.png")
plt.show()
In [62]:
loggs = np.arange(3.5, 5.6, 0.5)
In [63]:
flux_cont = [load_pixel(cont, 6000, logg, 0.0) for logg in loggs]
flux_mg_bot = [load_pixel(mg_bot, 6000, logg, 0.0) for logg in loggs]
flux_mg_side = [load_pixel(mg_side, 6000, logg, 0.0) for logg in loggs]
flux_fe_bot = [load_pixel(fe_bot, 6000, logg, 0.0) for logg in loggs]
flux_fe_side = [load_pixel(fe_side, 6000, logg, 0.0) for logg in loggs]
In [66]:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.plot(loggs, flux_cont, label="continuum")
ax.plot(loggs, flux_cont, "ko")
ax.plot(loggs, flux_mg_bot, label="Mg bottom")
ax.plot(loggs, flux_mg_bot, "ko")
ax.plot(loggs, flux_mg_side, label="Mg side")
ax.plot(loggs, flux_mg_side, "ko")
ax.plot(loggs, flux_fe_bot, label="Fe bottom")
ax.plot(loggs, flux_fe_bot, "ko")
ax.plot(loggs, flux_fe_side, label="Fe side")
ax.plot(loggs, flux_fe_side, "ko")
ax.set_xlabel(r"Temperature ($K$)")
ax.set_ylabel(r"$\propto f_\lambda$")
ax.legend(loc="lower right")
plt.show()
In [16]:
Zs = myInterface.points["Z"]
In [21]:
flux_cont = [load_pixel(cont, 6000, 4.0, Z) for Z in Zs]
flux_mg_bot = [load_pixel(mg_bot, 6000, 4.0, Z) for Z in Zs]
flux_mg_side = [load_pixel(mg_side, 6000, 4.0, Z) for Z in Zs]
flux_fe_bot = [load_pixel(fe_bot, 6000, 4.0, Z) for Z in Zs]
flux_fe_side = [load_pixel(fe_side, 6000, 4.0, Z) for Z in Zs]
In [23]:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.plot(Zs, flux_cont, label="continuum")
ax.plot(Zs, flux_cont, "ko")
ax.plot(Zs, flux_mg_bot, label="Mg bottom")
ax.plot(Zs, flux_mg_bot, "ko")
ax.plot(Zs, flux_mg_side, label="Mg side")
ax.plot(Zs, flux_mg_side, "ko")
ax.plot(Zs, flux_fe_bot, label="Fe bottom")
ax.plot(Zs, flux_fe_bot, "ko")
ax.plot(Zs, flux_fe_side, label="Fe side")
ax.plot(Zs, flux_fe_side, "ko")
ax.set_xlabel(r"Metallicity")
ax.set_ylabel(r"$\propto f_\lambda$")
ax.legend(loc="lower right")
plt.show()
In [22]:
flux_fe_bot
Out[22]:
In [13]:
from scipy.interpolate import RectBivariateSpline
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline as IUS
In [23]:
temps = myInterface.points["temp"]
print(temps)
loggs = myInterface.points["logg"]
print(loggs)
In [32]:
#This is the part that actually takes a while, loading all of the spectra.
#If we could do it on a computer with a lot of memory, it might not
#be so bad
values = np.empty((temps.size, loggs.size))
for i,temp in enumerate(temps):
for j, logg in enumerate(loggs):
values[i,j] = load_pixel(0, temp, logg, Z=0.0)
In [25]:
def interpolate_pixel(index, mytemp, mylogg):
#Create a 2D array of shape (temps.size, loggs.size) with the value of the flux at the pixel
#values = np.empty((temps.size, loggs.size))
#for i,temp in enumerate(temps):
# for j, logg in enumerate(loggs):
# values[i,j] = load_pixel(index, temp, logg, Z=0.0)
#create a bivariate spline over this grid
spl = RectBivariateSpline(temps, loggs, values)
return spl(mytemp, mylogg)
In [36]:
timeit value = interpolate_pixel(0, 6010, 3.8)
In [35]:
print(value)
In [12]:
plt.imshow(values,origin="upper",interpolation="none")
plt.show()
In [ ]: